Code
import dataclasses
import numpy as np
from typing import List, Tuple, UnionIn this tutorial, parameters of the model are stored in a class called NeuronParameters which is defined as follows:
@dataclasses.dataclass
class NeuronParameters:
C_m: float # membrane capacitance
E_L: float # leak reversal potential
E_K: float # potassium reversal potential
R_m: float # membrane resistance
V_th: float # threshold potential
V_reset: float # reset potential
V_peak: float = None # peak potential
# refractory period parameters
tau_vc: float = None # time constant for voltage clamp
tau_vth: float = None # time constant for voltage increase method
V_th_max: float = None # maximum voltage for V inc method
tau_G_ref: float = None # time constant for conductance increase method
# spike adaptation parameters
G_SRA: float = None # spike rate adaptation conductance
Delta_G_SRA: float = None # spike rate adaptation conductance delta
tau_SRA: float = None # spike rate adaptation time constant
V_max: float = None # maximum voltage for spike adaptation
Delta_th: float = None # spike adaptation threshold delta
a: float = None # spike adaptation parameter
b: float = None # spike adaptation parameter
def __post_init__(self):
self.G_L = 1./self.R_m # leak conductance
self.tau_m = self.C_m/self.G_L # membrane time constant
self.I_th = self.G_L*(self.V_th - self.E_L) # threshold current For simulation, we define a base class called SimulatorBase which is inherited by all the models. The base class defines the following attributes:
@dataclasses.dataclass
class SimulatorBase:
dt: float
times: np.ndarray
neuron_parameters: NeuronParameters
I: np.ndarray
num_fire: int = 0
def __post_init__(self):
self.V = np.zeros_like(self.times)
self.V[0] = self.neuron_parameters.E_L
@property
def fire_rate(self):
return self.num_fire/(self.times[-1] - self.times[0])
def dvdt(self):
raise NotImplementedError
def simulate(self):
raise NotImplementedErrorThe LIF model is given by the following equation:
C_m \frac{dV_m}{dt} = (E_L - V_m)/R_m + I_{app}
with the condition \quad V_m > V_{th} then V_m \rightarrow V_{reset}
C_m = 2.e-9 # membrane capacitance
E_L = -70.e-3 # leak reversal potential
E_K = -80e-3 # potassium reversal potential
R_m = 5.e6 # membrane resistance
G_L = 1./R_m # leak conductance
V_th = -50.e-3 # threshold potential
V_reset = -65.e-3 # reset potential
tau_m = C_m/G_L # membrane time constant
params = NeuronParameters(C_m=C_m, E_L=E_L, E_K=E_K, R_m=R_m, V_th=V_th, V_reset=V_reset)We can define a class called LIF which inherits from SimulatorBase and implements the dvdt and simulate methods:
@dataclasses.dataclass
class LIF(SimulatorBase):
noise_sigma: float = 0.
noises: np.ndarray = None
def __post_init__(self):
super().__post_init__()
self.noises = np.random.normal(size=self.times.shape) * self.noise_sigma * np.sqrt(self.dt)
def dvdt(self, V_m, I_app):
_dvdt = (self.neuron_parameters.E_L - V_m)/self.neuron_parameters.R_m + I_app
_dvdt = _dvdt/self.neuron_parameters.C_m
return _dvdt
def simulate(self):
for i in range(1, len(self.times)):
V_new = self.V[i-1] + self.dvdt(self.V[i-1], self.I[i-1]) * self.dt
if isinstance(self.noises, np.ndarray):
V_new = V_new + self.noises[i]
if V_new > self.neuron_parameters.V_th:
V_new = self.neuron_parameters.V_reset
self.num_fire += 1
self.last_fire_time = self.times[i]
self.V[i] = V_newWe can find the threshold current I_{th} by setting dV/dt of LIF equation to 0, giving the following equation: I_{th} = G_L(V_{th} - E_L)
Which is
In the plots below, the equation for the current threshold is validated; we can see that for 1% below the threshold current I_{th} no spikes are generated, but when increased to higher than the threshold current, spikes can be seen.
import matplotlib.pyplot as plt
import numpy as np
import plotly.express as px
from plotly.subplots import make_subplots
import plotly.graph_objects as go
dt = 0.1e-4
times = np.arange(0, 200.e-3, dt)
LIF_below_th = LIF(dt, times, params, np.ones_like(times)*0.99*params.I_th)
LIF_below_th.simulate()
LIF_above_th = LIF(dt, times, params, np.ones_like(times)*1.01*params.I_th)
LIF_above_th.simulate()
# plot
fig = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=times, y=LIF_below_th.V, name='0.99 * I_th'), row=1, col=1)
fig.add_trace(go.Scatter(x=times, y=LIF_above_th.V, name='1.01 * I_th'), row=2, col=1)
fig.update_layout(height=600, width=800, title_text="Voltage vs time")
# set axis names
fig.update_xaxes(title_text="time (s)", row=2, col=1)
fig.update_yaxes(title_text="Voltage (V)", row=1, col=1)
fig.update_yaxes(title_text="Voltage (V)", row=2, col=1)
fig.show()The firing rate of the neuron is given by the following equation:
1/f = \tau_m \ln \left(\frac{I R_m + E_L - V_{reset}}{I R_m + E_L - V_{th}}\right)
Below plotting the equation for the firing rate of the neuron using the above equation and the firing rate of the neuron using the LIF equation. We can see that the firing rate of the neuron is very close to the firing rate of the neuron using the LIF equation. We can also see that the firing rate is 0 when the current is below the threshold current I_{th}=4nA.
fr_array = []
fr_eq_array = []
I_array = np.arange(2e-9, 1e-8, 1e-9)
def fr_func(I):
# The equation for the firing rate of the neuron
value = params.tau_m * np.log(I * params.R_m + params.E_L - params.V_reset) - params.tau_m * np.log(I*params.R_m + params.E_L-params.V_th)
if value < 0:
return 0
value = 1/value
return value
for I in I_array:
with np.errstate(divide='ignore', invalid='ignore'):
# from sim
times = np.arange(0, 200.e-3, dt)
LIF_model = LIF(dt, times, params, np.ones_like(times)*I)
LIF_model.simulate()
fr_array.append(LIF_model.fire_rate)
# From eq
fr_eq = fr_func(I)
fr_eq_array.append(fr_eq)
# use plotly
fig = make_subplots(rows=1, cols=1, shared_xaxes=True, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=I_array, y=fr_array, name='f LIF'), row=1, col=1)
fig.add_trace(go.Scatter(x=I_array, y=fr_eq_array, name='f equation'), row=1, col=1)
fig.update_layout(height=600, width=800, title_text="f vs I")
# set axis names
fig.update_xaxes(title_text="I (A)", row=1, col=1)
fig.update_yaxes(title_text="f (Hz)", row=1, col=1)
fig.show()Using the Euler-Mayamara method, we add random noise to the current by \sigma w(t); where w(t) is a white noise generator with zero mean and unit variance, and \sigma scales the noise.
dt = 0.1e-3
times = np.arange(0, 150.e-3, dt)
fr_array_no_noise = []
fr_array_noise1 = []
fr_array_noise2 = []
I_array = np.arange(1e-9, 9e-9, 1e-10)
sigma_1 = 1e-2
sigma_2 = 3e-2
for I in I_array:
LIF_no_noise = LIF(dt, times, params, np.ones_like(times)*I)
LIF_noise1 = LIF(dt, times, params, np.ones_like(times)*I, noise_sigma=sigma_1)
LIF_noise2 = LIF(dt, times, params, np.ones_like(times)*I, noise_sigma=sigma_2)
LIF_no_noise.simulate()
LIF_noise1.simulate()
LIF_noise2.simulate()
fr_array_no_noise.append(LIF_no_noise.fire_rate)
fr_array_noise1.append(LIF_noise1.fire_rate)
fr_array_noise2.append(LIF_noise2.fire_rate)
# use plotly
fig = make_subplots(rows=1, cols=1, shared_xaxes=True, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=I_array, y=fr_array_no_noise, name='sigma_I = 0'), row=1, col=1)
fig.add_trace(go.Scatter(x=I_array, y=fr_array_noise1, name=f'sigma_I={sigma_1}'), row=1, col=1)
fig.add_trace(go.Scatter(x=I_array, y=fr_array_noise2, name=f'sigma_I={sigma_2}'), row=1, col=1)
fig.update_layout(height=600, width=800, title_text="f vs I")
# set axis names
fig.update_xaxes(title_text="I (A)", row=1, col=1)
fig.update_yaxes(title_text="f (Hz)", row=1, col=1)
fig.show()In the above figure \sigma_I takes the values [0, 1e^{-2}, 3e^{-2}]. It seems that by increasing the noise, the curve is smoothed out near the threshold current region. We can see that the blue curve (no noise) has a near verticle section near the threshold current, but as we increase the noise, the curve is more “smooth” around the threshold current I_{th}. I suppose this is because in real-life situations the rather sharp threshold current is not exactly the same as the threshold current in the model. This is because the model is a simplification of the real-life situation.
I understand the basic properties of the leaky integrate-and-fire (LIF) neuron model and how noise causes the smoothing of f-I curves.
In this tutorial, I looked at 3 different ways to implement the refractory period in the LIF neuron model.
@dataclasses.dataclass
class LIF_Voltage_Clamp(LIF):
# LIF with voltage clamp
last_spike_time: float = -np.inf # Last spike time
spike_locations: List[int] = dataclasses.field(default_factory=list) # Spike locations
def simulate(self):
for i in range(1, len(self.times)):
t = self.times[i]
if t - self.last_spike_time < self.neuron_parameters.tau_vc:
self.V[i] = self.neuron_parameters.V_reset
continue
V_new = self.V[i-1] + self.dvdt(self.V[i-1], self.I[i-1]) * self.dt
V_new = V_new + self.noises[i]
if V_new > self.neuron_parameters.V_th:
self.num_fire += 1
self.V[i] = self.neuron_parameters.V_reset
self.last_spike_time = t
self.spike_locations.append(i)
else:
self.V[i] = V_new
# set V_peak
self.V[self.spike_locations] = self.neuron_parameters.V_peakThis method is a combination of the voltage clamp method and the threshold increase method. The refractory conductance is added to the voltage clamp method.
@dataclasses.dataclass
class LIF_Threshold_Increase(LIF):
use_refractory_conductance: bool = False # Use refractory conductance
V_th_array: np.ndarray = None # Threshold voltage array
spike_locations: List[int] = dataclasses.field(default_factory=list) # Spike locations
# refractory conductance
G_ref_array: np.ndarray = None # Refractory conductance array
def __post_init__(self):
super().__post_init__()
self.V_th_array = np.ones_like(self.times) * self.neuron_parameters.V_th
self.V_th_array[0] = self.neuron_parameters.V_th
self.G_ref_array = np.zeros_like(self.times)
def dv_thdt(self, V_th):
_dv_thdt = (self.V_th_array[0] - V_th) / self.neuron_parameters.tau_vth
return _dv_thdt
def dG_refdt(self, G_ref):
_dG_refdt = - G_ref / self.neuron_parameters.tau_G_ref
return _dG_refdt
def dvdt(self, V, I, time_index):
_dvdt = super().dvdt(V, I)
if self.use_refractory_conductance:
# add refractory conductance
_dvdt = _dvdt + self.G_ref_array[time_index] * ( self.neuron_parameters.E_K - V) / self.neuron_parameters.C_m
self.G_ref_array[time_index+1] = self.G_ref_array[time_index] + self.dG_refdt(self.G_ref_array[time_index]) * self.dt
return _dvdt
def simulate(self):
for i in range(1, len(self.times)):
time = self.times[i]
V_new = self.V[i-1] + self.dvdt(self.V[i-1], self.I[i-1], i-1) * self.dt
V_new = V_new + self.noises[i]
self.V[i] = V_new
self.V_th_array[i] = self.V_th_array[i-1] + self.dv_thdt(self.V_th_array[i-1]) * self.dt
if V_new > self.V_th_array[i]:
self.spike_locations.append(i)
self.V_th_array[i] = self.neuron_parameters.V_th_max
self.num_fire += 1
if self.use_refractory_conductance:
# add 2 micro Siemens to the refractory conductance
self.G_ref_array[i] = self.G_ref_array[i] + 2.e-6
else:
self.V[i] = self.neuron_parameters.V_reset
#set V_peak
self.V[self.spike_locations] = self.neuron_parameters.V_peakdt = 1e-5
times = np.arange(0, 100.e-3, dt)
currents_1 = np.ones_like(times)*220e-12
lif_vc1 = LIF_Voltage_Clamp(dt, times, params, currents_1, noise_sigma=1.5e-2)
lif_vc1.simulate()
currents_2 = np.ones_like(times)*600e-12
lif_vc2 = LIF_Voltage_Clamp(dt, times, params, currents_2, noise_sigma=1.5e-2)
lif_vc2.simulate()
# use plotly to plot the two simulations with vc on the same trace
fig = make_subplots(rows=1, cols=1, shared_xaxes=True, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=times, y=lif_vc1.V, name='I=220pA'), row=1, col=1)
fig.add_trace(go.Scatter(x=times, y=lif_vc2.V, name='I=600pA'), row=1, col=1)
fig.update_layout(height=600, width=800, title_text="Forced VC: Voltage vs Time")
# set axis names
fig.update_xaxes(title_text="Time (s)", row=1, col=1)
fig.update_yaxes(title_text="Voltage (V)", row=1, col=1)
fig.show()lif_ts1 = LIF_Threshold_Increase(dt, times, params, currents_1, noise_sigma=1.5e-2)
lif_ts1.simulate()
lif_ts2 = LIF_Threshold_Increase(dt, times, params, currents_2, noise_sigma=1.5e-2)
lif_ts2.simulate()
# use plotly to plot the two simulations on the same trace
fig = make_subplots(rows=1, cols=1, shared_xaxes=True, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=times, y=lif_ts1.V, name='I=220pA'), row=1, col=1)
fig.add_trace(go.Scatter(x=times, y=lif_ts2.V, name='I=600pA'), row=1, col=1)
fig.add_trace(go.Scatter(x=times, y=lif_ts1.V_th_array, name='V_th for I=220pA', line=dict(color='green', dash='dash')), row=1, col=1)
fig.add_trace(go.Scatter(x=times, y=lif_ts2.V_th_array, name='V_th for I=600pA', line=dict(color='cyan', dash='dash')), row=1, col=1)
fig.update_layout(height=600, width=800, title_text="Threshold increase: Voltage vs Time")
# set axis names
fig.update_xaxes(title_text="Time (s)", row=1, col=1)
fig.update_yaxes(title_text="Voltage (V)", row=1, col=1)
fig.show()lif_rc1 = LIF_Threshold_Increase(dt, times, params, currents_1, noise_sigma=1.5e-2, use_refractory_conductance=True)
lif_rc1.simulate()
lif_rc2 = LIF_Threshold_Increase(dt, times, params, currents_2, noise_sigma=1.5e-2, use_refractory_conductance=True)
lif_rc2.simulate()
# use plotly to plot the two simulations on the same trace
fig = make_subplots(rows=1, cols=1, shared_xaxes=True, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=times, y=lif_rc1.V, name='I=220pA'), row=1, col=1)
fig.add_trace(go.Scatter(x=times, y=lif_rc2.V, name='I=600pA'), row=1, col=1)
fig.update_layout(height=600, width=800, title_text="Threshold increase with refractory conductance: Voltage vs Time")
# set axis names
fig.update_xaxes(title_text="Time (s)", row=1, col=1)
fig.update_yaxes(title_text="Voltage (V)", row=1, col=1)
fig.show()currents = np.arange(100e-12, 600e-12, 20e-12)
dt = 1e-4
times = np.arange(0., 0.2, dt)
# create empty lists to store the firing rates
fr_fc = []
fr_ts = []
fr_rcts = []
# mean V vs different currents
V_fc = []
V_ts = []
V_rcts = []
for current in currents:
lif_fc = LIF_Voltage_Clamp(dt, times, params, np.ones_like(times)*current, noise_sigma=1.e-3)
lif_fc.simulate()
lif_ts = LIF_Threshold_Increase(dt, times, params, np.ones_like(times)*current, noise_sigma=1.e-3)
lif_ts.simulate()
lif_rcts = LIF_Threshold_Increase(dt, times, params, np.ones_like(times)*current, noise_sigma=1.e-3, use_refractory_conductance=True)
lif_rcts.simulate()
fr_fc.append(lif_fc.fire_rate)
fr_ts.append(lif_ts.fire_rate)
fr_rcts.append(lif_rcts.fire_rate)
V_fc.append(np.mean(lif_fc.V))
V_ts.append(np.mean(lif_ts.V))
V_rcts.append(np.mean(lif_rcts.V))The following plots are shown:
# Plot the three refractory models on the same plot, don't use shared x axes and legend
fig = make_subplots(rows=1, cols=1, shared_xaxes=False, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=currents, y=fr_fc, name='Voltage clamp'), row=1, col=1)
fig.add_trace(go.Scatter(x=currents, y=fr_ts, name='Threshold increase'), row=1, col=1)
fig.add_trace(go.Scatter(x=currents, y=fr_rcts, name='Threshold increase with refractory conductance'), row=1, col=1)
fig.update_layout(height=600, width=800, title_text="Mean firing rate vs input current")
# set axis names
fig.update_xaxes(title_text="Input current (A)", row=1, col=1)
fig.update_yaxes(title_text="Mean firing rate (Hz)", row=1, col=1)
fig.show()# Plot the three refractory models on the same plot, don't use shared x axes and legend
fig = make_subplots(rows=1, cols=1, shared_xaxes=False, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=currents, y=V_fc, name='Voltage clamp'), row=1, col=1)
fig.add_trace(go.Scatter(x=currents, y=V_ts, name='Threshold increase'), row=1, col=1)
fig.add_trace(go.Scatter(x=currents, y=V_rcts, name='Threshold increase with refractory conductance'), row=1, col=1)
fig.update_layout(height=600, width=800, title_text="Mean voltage vs input current")
# set axis names
fig.update_xaxes(title_text="Input current (A)", row=1, col=1)
fig.update_yaxes(title_text="Mean voltage (V)", row=1, col=1)
fig.show()Upon examining the mean voltage versus input current plot shown above, we can observe a non-monotonic relationship when the threshold increase method is employed. In the case of threshold increase, the voltage threshold decreases at the same rate as the input current increases, causing the voltage to reach the threshold more quickly and higher than V_{th}. This explains why the mean voltage is not monotonically increasing in the threshold increase case. As we increase the input current, the mean voltage decreases for a while, then increases again.
However, when we add refractory conductance to the threshold increase method, the mean voltage decreases monotonically as the input current increases. This suggests that the refractory conductance may be slowing down the voltage from reaching the threshold.
# Plot the three refractory models on the same plot, don't use shared x axes and legend
fig = make_subplots(rows=1, cols=1, shared_xaxes=False, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=fr_fc, y=V_fc, name='Voltage clamp'), row=1, col=1)
fig.add_trace(go.Scatter(x=fr_ts, y=V_ts, name='Threshold increase'), row=1, col=1)
fig.add_trace(go.Scatter(x=fr_rcts, y=V_rcts, name='Threshold increase with refractory conductance'), row=1, col=1)
fig.update_layout(height=600, width=800, title_text="Mean voltage vs mean firing rate")
# set axis names
fig.update_xaxes(title_text="Mean firing rate (Hz)", row=1, col=1)
fig.update_yaxes(title_text="Mean voltage (V)", row=1, col=1)
fig.show()In this tutorial, I learnt about the different ways to model the refractory period using 3 methods. The non-monotonic relationship for the mean voltage versus input current plot was interesting to observe. My hypothesis was that the
The adaptionation currents adds to the LIF model by adding a hyperpolarising current.
C_m \frac{dV}{dt} = (E_L - V) / R_m + G_{SRA} (E_{SRA} - V) + I_{app}
where
\frac{G_{SRA}}{dt} = -G_{SRA} / \tau_{SRA}
and,
\quad V \ge V_{th} \implies V \rightarrow V_{TH} \And G_{SRA} \rightarrow G_{SRA} + \Delta G_{SRA}
@dataclasses.dataclass
class LIF_SRA(LIF):
# LIF model with spike-rate adaptation
G_SRA_array: np.ndarray = None # spike-rate adaptation conductance
last_spike_time: float = None # time of last spike
isi_array: List[float] = dataclasses.field(default_factory=list) # inter-spike interval
def __post_init__(self):
super().__post_init__()
self.G_SRA_array = np.zeros(self.times.shape)
def dvdt(self, V_m, I_app, time_index):
lif_dvdt = super().dvdt(V_m, I_app)
# spike-rate adaptation
G_SRA = self.G_SRA_array[time_index]
dG_SRA = -G_SRA / self.neuron_parameters.tau_SRA
self.G_SRA_array[time_index + 1] = G_SRA + dG_SRA * self.dt
return lif_dvdt + G_SRA * (self.neuron_parameters.E_K - V_m)/self.neuron_parameters.C_m
def simulate(self):
for i in range(1, len(self.times)):
self.V[i] = self.V[i - 1] + self.dvdt(self.V[i - 1], self.I[i - 1], i - 1) * self.dt
if self.V[i] >= self.neuron_parameters.V_th:
self.V[i] = self.neuron_parameters.V_reset
self.G_SRA_array[i] += self.neuron_parameters.Delta_G_SRA
self.num_fire += 1
if self.last_spike_time is not None:
self.isi_array.append(self.times[i] - self.last_spike_time)
self.last_spike_time = self.times[i]Below we plot the current, voltage and spike-rate adaptation conductance for a single neuron with a constant input current of 500 pA.
params = NeuronParameters(
C_m=100e-12,
E_L=-75.e-3,
R_m=100.e6,
E_K=-80.e-3,
V_th=-50e-3,
V_reset=-80.e-3,
tau_vc=5.e-3,
V_peak=50.e-3,
tau_SRA=200e-3,
Delta_G_SRA=1e-9,
)
dt = 1e-3
times = np.arange(0, 1.5, dt)
currents = np.zeros_like(times)
mask = np.logical_and(times > 0.5, times < 1.0)
currents[mask] = 400e-12
lif_sra = LIF_SRA(dt, times, params, currents, noise_sigma=0.0)
lif_sra.simulate()
fig = make_subplots(rows=1, cols=1, shared_xaxes=False, vertical_spacing=0.02)
fig.add_trace(go.Scatter(x=lif_sra.times, y=lif_sra.V, name='Voltage'), row=1, col=1)
fig.add_trace(go.Scatter(x=lif_sra.times, y=lif_sra.G_SRA_array, name='G_SRA(t)'), row=1, col=1)
# Plot voltage and current, voltage and spike-rate adaptation conductance in different plots
# 3 rows using subplots
fig = make_subplots(rows=4, cols=1, shared_xaxes=False, vertical_spacing=0.02)
# plot current
fig.add_trace(go.Scatter(x=lif_sra.times, y=lif_sra.I, name='Current'), row=1, col=1)
# plot voltage
fig.add_trace(go.Scatter(x=lif_sra.times, y=lif_sra.V, name='Voltage'), row=2, col=1)
# plot spike-rate adaptation conductance
fig.add_trace(go.Scatter(x=lif_sra.times, y=lif_sra.G_SRA_array, name='G_SRA(t)'), row=3, col=1)
# plot the IST
#fig.add_trace(go.Scatter(x=lif_sra.times, y=lif_sra.inter_spike_intervals, name='ISI(t)'), row=4, col=1)
# set axis names
fig.update_xaxes(title_text="Time (s)", row=3, col=1)
fig.update_yaxes(title_text="Current (A)", row=1, col=1)
fig.update_yaxes(title_text="Voltage (V)", row=2, col=1)
#fig.update_yaxes(title_text="G_SRA (S)", row=3, col=1)
fig.show()Wow, the above plot is really interesting and pretty. The adaptation current - an outward hyperpolarising current - increases when spikes are fired and decrease after due to the differential equation, but when another spike occurs, the adaptation current increases again. This adaptation current is a function of the membrane potential meaning that the neuron is adapting to the firing rate of the neuron.
currents = np.arange(150e-12, 350e-12, 10e-12)
times = np.arange(0, 5, dt)
first_isi = []
steady_isi = []
single_spike = []
for i, I in enumerate(currents):
lif_sra = LIF_SRA(dt, times, params, I * np.ones_like(times), noise_sigma=1e-5)
lif_sra.simulate()
isi_array = lif_sra.isi_array
if lif_sra.num_fire > 3:
first_isi.append(1 / np.mean(isi_array[0:3]))
steady_isi.append(1 / isi_array[-1])
elif lif_sra.num_fire == 1:
single_spike.append(i)
else:
first_isi.append(0)
steady_isi.append(0)
fig = go.Figure()
# plot first with circles
fig.add_trace(go.Scatter(x=currents, y=first_isi, mode='markers', name='1/First ISI'))
fig.add_trace(go.Scatter(x=currents, y=steady_isi, name='1/Steady ISI'))
# plot single spike with X
fig.add_trace(go.Scatter(x=currents[single_spike], y=np.zeros_like(single_spike), mode='markers', name='Single spike'))
fig.update_xaxes(title_text="Current (A)")
fig.update_yaxes(title_text="Firing rate (Hz)")
fig.show()In the above f-i curve, the inverse first ISI and inverse steady state ISI shows similarity near the threshold current, and as the applied current increases the inverse first ISI increases faster than inv steady state ISI. This is because due to the spike-rate adaptation, the neuron initially fires at a higher rate, but at steady state the neuron fires much slower.
@dataclasses.dataclass
class AELIF(LIF):
I_SRA_array: np.ndarray = None
last_spike_time: float = None # time of last spike
isi_array: List[float] = dataclasses.field(default_factory=list) # inter-spike interval
def __post_init__(self):
super().__post_init__()
self.I_SRA_array = np.zeros_like(self.times)
def dI_SRA_dt(self, V_m, I_SRA):
a = self.neuron_parameters.a
b = self.neuron_parameters.b
_dI_SRA_dt = a * (V_m - self.neuron_parameters.E_L) - I_SRA
return _dI_SRA_dt/self.neuron_parameters.tau_SRA
def dvdt(self, V_m, I_app, I_SRA):
exp_term = self.neuron_parameters.Delta_th * np.exp((V_m - self.neuron_parameters.V_th)/self.neuron_parameters.Delta_th)
_dv_dt = self.neuron_parameters.G_L*(self.neuron_parameters.E_L - V_m + exp_term) - I_SRA + I_app
#print(f'V_m: {V_m}, I_app: {I_app}, I_SRA: {I_SRA}, dv_dt: {_dv_dt/self.neuron_parameters.C_m}')
return _dv_dt/self.neuron_parameters.C_m
def simulate(self):
for i in range(1, len(self.times)):
self.V[i] = self.V[i-1] + self.dt * self.dvdt(self.V[i-1], self.I[i-1], self.I_SRA_array[i-1])
self.I_SRA_array[i] = self.I_SRA_array[i-1] + self.dt * self.dI_SRA_dt(self.V[i-1], self.I_SRA_array[i-1])
if self.V[i] >= self.neuron_parameters.V_max:
self.V[i] = self.neuron_parameters.V_reset
self.I_SRA_array[i] += self.neuron_parameters.b
self.num_fire += 1
if self.last_spike_time is not None:
self.isi_array.append(self.times[i] - self.last_spike_time)
self.last_spike_time = self.times[i]Adding an exponential term to the adaptation current model we get the AELIF model:
C_m \frac{dV_m}{dt} = G_L(E_L - V_m + \Delta_{th}e^{\frac{V_m - V_{th}}{\Delta_{th}}}) - I_{SRA} + I_{app}
The exponential term creates a positive feedback loop, which makes the neuron membrane potential spike faster and faster. Using this we don’t need to have the IF statement when V_m > V_{th} \rightarrow V_{peak}, however this still has an IF statement because of the clamping by V_{max}, how can we get rid of this clamping?
params.a = 2e-9
params.b = 0.02e-9
params.Delta_th = 2e-3
params.tau_SRA = 200e-3
params.G_L = 10e-9
params.V_max = 100e-3
# plot for 500 pA
dt = 1e-3
times = np.arange(0, 1.5, dt)
currents = np.zeros_like(times)
mask = np.logical_and(times > 0.5, times < 1.0)
currents[mask] = 300e-12
lif_sra = AELIF(dt, times, params, currents, noise_sigma=0)
lif_sra.simulate()
# plot on different subplots
fig = make_subplots(rows=3, cols=1, subplot_titles=("AELIF model", "I_SRA"))
fig.add_trace(go.Scatter(x=times, y=lif_sra.V), row=1, col=1)
fig.add_trace(go.Scatter(x=times, y=lif_sra.I_SRA_array), row=2, col=1)
fig.add_trace(go.Scatter(x=times, y=currents), row=3, col=1)
fig.update_xaxes(title_text="Time (s)", row=3, col=1)
fig.update_yaxes(title_text="V_m (V)", row=1, col=1)
fig.update_yaxes(title_text="I_SRA (A)", row=2, col=1)
fig.update_yaxes(title_text="I_app (A)", row=3, col=1)
fig.update_layout(showlegend=False)
fig.show()
dt = 1e-4
currents = np.arange(150e-12, 350e-12, 5e-12)
times = np.arange(0, 5, dt)
first_isi = []
steady_isi = []
single_spike = []
for i, I in enumerate(currents):
lif_sra = AELIF(dt, times, params, I * np.ones_like(times), noise_sigma=0) #=1e-5)
lif_sra.simulate()
isi_array = lif_sra.isi_array
if lif_sra.num_fire > 3:
first_isi.append(1 / np.mean(isi_array[0:3]))
steady_isi.append(1 / isi_array[-1])
elif lif_sra.num_fire == 1:
single_spike.append(i)
else:
first_isi.append(0)
steady_isi.append(0)
print(single_spike)
fig = go.Figure()
# plot first with circles
fig.add_trace(go.Scatter(x=currents, y=first_isi, mode='markers', name='1/First ISI'))
fig.add_trace(go.Scatter(x=currents, y=steady_isi, name='1/Steady ISI'))
# plot single spike with X
fig.add_trace(go.Scatter(x=currents[single_spike], y=np.zeros_like(single_spike), mode='markers', name='Single spike'))
fig.update_xaxes(title_text="Current (A)")
fig.update_yaxes(title_text="Firing rate (Hz)")
fig.show()[19, 20, 21, 22]
Looking at the two f-I curves they seem to be quite similar, the AELIF is slightly shifted to the right and also the inverse of the first couple ISI is slightly lower. This might be because in the AELIF model the exponential term is a positive feed-back loop, but not instantaneous, so it takes a little bit of time for the neuron to spike, thus reducing the firing rate for the first couple of spikes.
In chapter 2 I learnt the following: